if !d.name eq 'PS' then begin
    device,xsize=20,ysize=9,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
@rhoprof
;!p.charsize=2.5
!p.multi=[0,4,1] & !p.charsize=2.0 
!x.margin=[7,0] & !y.margin=[5,3]

default,nn,5
ntimes = (size(u))[3]
good = (ntimes-nn)+indgen(nn)-1
print,'doing averages over',nn, ' xaver records ...'

ruu = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
rvv = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
rww = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)


uu=fltarr(l,m)
vv=fltarr(l,m)
ww=fltarr(l,m)
oox=fltarr(l,m)
ooy=fltarr(l,m)
ooz=fltarr(l,m)

omega=fltarr(l,m)

quu=fltarr(l,m)
qvv=fltarr(l,m)
qww=fltarr(l,m)
qwv=fltarr(l,m)
qwu=fltarr(l,m)
qvu=fltarr(l,m)
rsinth=fltarr(l,m)
rurms2=fltarr(l,m)
frmc=fltarr(l,m)
frrs=fltarr(l,m)
fthmc=fltarr(l,m)
fthrs=fltarr(l,m)

frt=fltarr(l,m)
ftht=fltarr(l,m)
flrs=fltarr(l,m)
flmc=fltarr(l,m)
fzrs=fltarr(l,m)
fzmc=fltarr(l,m)
fx=fltarr(l,m)
fy=fltarr(l,m)

urms2t = ruu+rvv+rww
urms2 = (total(urms2t[good]/float((size(good))[1])))
omega = fltarr(l,m)
omega_0=2.*!pi/(14.*24.*3600.) 
; reynolds stresses
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=rr*r[k]*sin(theta[j])
    omega[k,j] = omega_0*rsinth[k,j]
    uu[k,j] = total(u[j,k,good])/float((size(good))[1])
    vv[k,j] = total(v[j,k,good])/float((size(good))[1])
    ww[k,j] = total(w[j,k,good])/float((size(good))[1])
    oox[k,j] = total(ox[j,k,good])/float((size(good))[1])
    ooy[k,j] = total(oy[j,k,good])/float((size(good))[1])
    ooz[k,j] = total(oz[j,k,good])/float((size(good))[1])
    quu[k,j]=rho[j,k]*total(u2[j,k,good]- u[j,k,good]^2)/float((size(good))[1])
    qvv[k,j]=rho[j,k]*total(v2[j,k,good]- v[j,k,good]^2)/float((size(good))[1])
    qww[k,j]=rho[j,k]*total(w2[j,k,good]- w[j,k,good]^2)/float((size(good))[1])
;
    qwv[k,j]=total(rwv[j,k,good]- rho[j,k]*oz[j,k,good]*v[j,k,good])/float((size(good))[1])
    qwu[k,j]=total(rwu[j,k,good]- rho[j,k]*oz[j,k,good]*u[j,k,good])/float((size(good))[1])
    qvu[k,j]=total(rvu[j,k,good]- rho[j,k]*oy[j,k,good]*u[j,k,good])/float((size(good))[1])
;   
    qwv[k,j]=rh[k]*qwv[k,j]/rho[j,k]
    qwu[k,j]=rh[k]*qwu[k,j]/rho[j,k]
    qvu[k,j]=rh[k]*qvu[k,j]/rho[j,k]
;
    rurms2[k,j]=quu[k,j]+qvv[k,j]+qww[k,j]
; r components of the angular momentum flux. all normalized with urms2
    frmc[k,j]=rsinth[k,j]*rh[k]*(uu[k,j])*ww[k,j]
    frrs[k,j]=rsinth[k,j]*qwu[k,j]
; th components of the angular momentum flux. all normalized with urms2
    fthmc[k,j] = rsinth[k,j]*rh[k]*(uu[k,j])*vv[k,j]
    fthrs[k,j] = rsinth[k,j]*qvu[k,j]
; totals    
    frt[k,j] =  (frmc[k,j]+frrs[k,j])
    ftht[k,j] = (fthmc[k,j]+fthrs[k,j])
;
; angular momentum flux in cilindrical radius \lambda
    flrs[k,j]=frrs[k,j]*sin(theta[j]) + fthrs[k,j]*cos(theta[j])
    flmc[k,j]=frmc[k,j]*sin(theta[j]) + fthmc[k,j]*cos(theta[j])
;
    fzrs[k,j]=frrs[k,j]*cos(theta[j]) - fthrs[k,j]*sin(theta[j])
    fzmc[k,j]=frmc[k,j]*cos(theta[j]) - fthmc[k,j]*sin(theta[j])

;
;   in order to make arrow plots
;    fx[k,j]=frt[k,j]*sin(theta[j]) + ftht[k,j]*cos(theta[j])
;    fy[k,j]=frt[k,j]*cos(theta[j]) - ftht[k,j]*sin(theta[j])
;
;; terms to compute radial derivative
;    r2frmc[k,j]=(rr*r[k])^2*frmc[k,j]
;    r2frrs[k,j]=(rr*r[k])^2*frrs[k,j]
;; terms to compute latitudinal derivative
;    sfthmc[k,j]=sin(theta[j]*fthmc[k,j]
;    sfthrs[k,j]=sin(theta[j]*fthrs[k,j]
  endfor
endfor

; radial derivatives
;for j=0, m-1 do begin
;  divrMC[*,j]=rr*deriv(r[*],r2frmc[*,j])
;  divrRS[*,j]=rr*deriv(r[*],r2frrs[*,j])
;endfor
;for k=0, l-1 do begin
;  divthMC[k,*]=
;endfor

;varpi_letter="166B
;varpi='!9'+String(varpi_letter)+'!X'
varpi='!S!7x!R!A_!N!X'

flrs=flrs/1e12
flmc=flmc/1e12
fzrs=fzrs/1e12
fzmc=fzmc/1e12

xlength=0.19
xinterval=0.05
bini=0.10
bint=0.24
xpos1=xinterval
xpos2=xinterval+xlength

cpos=[xpos1,0.04,xpos2,0.94]
bpos=[bini,0.31,bini+0.01,0.63]
lev = 2.*max(abs(flrs))*indgen(nlev)/float(nlev-1)-max(abs(flrs))
;lev=-1.1+indgen(64)*1.6/63
contour,smooth(flrs[*,thg],0),x[*,thg],y[*,thg],/fi,nl=64,/iso,pos=cpos,$
;	title='!8F'+varpi+'!N!6 ',$
        lev=lev,xtitle='!8x!6',ytitle='!8y!6',xstyle=4,ystyle=4
colorbar,range=[min(lev),max(lev)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.7,title='!6RS (!810!U12!N !6kg s!U-2!N!6) '
plots,circle(0,0,0.96),thick=6
plots,circle(0,0,0.61),thick=6
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=6,li=0
oplot,[0,0],[0.61,0.96],thick=6,li=0
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
bini=bini+bint

lev = 2.*max(abs(flmc))*indgen(nlev)/float(nlev-1)-max(abs(flmc))
cpos=[xpos1,0.04,xpos2,0.94]
bpos=[bini,0.31,bini+0.01,0.63]
jj=thg
contour,smooth(flmc[*,jj],2),x[*,jj],y[*,jj],/fi,nl=nlev,/iso,pos=cpos,$
;        title='!8F'+varpi+'!6 ',$
	xtitle='!8x!6',xstyle=4,ystyle=4,lev=lev
colorbar,range=[min(lev),max(lev)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.7,title='!6MC (!810!U12!N !6kg s!U-2!N!6) '
jj=where(th GE -88 AND th LT 88)
;partvelvec,smooth(ux[*,jj],2),smooth(uy[*,jj],2),x[*,jj],y[*,jj],/over,fraction=0.08,length=0.1
plots,circle(0,0,0.96),thick=6
plots,circle(0,0,0.61),thick=6
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=6,li=0
oplot,[0,0],[0.61,0.96],thick=6,li=0
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
bini=bini+bint
print,bini

lev = 2.*max(abs(fzrs))*indgen(nlev)/float(nlev-1)-max(abs(fzrs))
cpos=[xpos1,0.04,xpos2,0.94]
bpos=[bini,0.31,bini+0.01,0.63]
jj=thg
contour,smooth(fzrs[*,jj],2),x[*,jj],y[*,jj],/fi,nl=nlev,/iso,pos=cpos,$
;        title='!8Fz!N!6 ',$
	xtitle='!8x!6',xstyle=4,ystyle=4,lev=levi
colorbar,range=[min(lev),max(lev)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.7,title='!6RS (!810!U12!N !6kg s!U-2!N!6) '
jj=where(th GE -88 AND th LT 88)
;partvelvec,smooth(ux[*,jj],2),smooth(uy[*,jj],2),x[*,jj],y[*,jj],/over,fraction=0.08,length=0.1
plots,circle(0,0,0.96),thick=6
plots,circle(0,0,0.61),thick=6
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=6,li=0
oplot,[0,0],[0.61,0.96],thick=6,li=0
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
bini=bini+bint

lev = 2.*max(abs(fzmc))*indgen(nlev)/float(nlev-1)-max(abs(fzmc))
cpos=[xpos1,0.04,xpos2,0.94]
bpos=[bini,0.31,bini+0.01,0.63]
jj=thg
contour,smooth(fzmc[*,jj],2),x[*,jj],y[*,jj],/fi,nl=nlev,/iso,pos=cpos,$
;        title='!8Fz!N!6 ',$
	xtitle='!8x!6',xstyle=4,ystyle=4,lev=lev
colorbar,range=[min(lev),max(lev)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.7,title='!6MC (!810!U12!N !6kg s!U-2!N!6) '
jj=where(th GE -88 AND th LT 88)
;partvelvec,smooth(ux[*,jj],2),smooth(uy[*,jj],2),x[*,jj],y[*,jj],/over,fraction=0.08,length=0.1
plots,circle(0,0,0.96),thick=6
plots,circle(0,0,0.61),thick=6
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=6,li=0
oplot,[0,0],[0.61,0.96],thick=6,li=0


!p.multi=0
end
